library(ggplot2)
# Set working directory to project directory
setwd(wd)

load("data/dataforstage2.Rda")

# Issue voting over time
# Identify five issues with highest scores
#  Figure 1a
df <- dat[order(dat$iv_changes_gov_abs,decreasing=T),]
df$issueshort <- NA
df$issueshort[1:5] <- c("Strengthen WEF", "Nationalization","Cancel WEF","Cancel WEF","Occupational Pension") 
pdf("absovertime.pdf")
ggplot(df,aes(year,iv_changes_gov_abs)) + geom_point() + theme_bw() + geom_label(data=subset(df,!is.na(df$issueshort)),mapping=aes(year,iv_changes_gov_abs + c(1.1,1.1,1.1,-1.1,1.1),label=issueshort)) + xlab("Year") + ylab("Absolute Issue Voting") + xlim(1954,2010)
dev.off()


# Individual data points. Figure 1b
pdf("ivpublicsupport.pdf")
ggplot(dat,aes(gw_fvr_dk,iv_changes_gov)) + geom_point() + theme_bw() + geom_hline(yintercept=0,linetype="dashed") + ylab("Total Issue Voting") + xlab("Public Support (%)")
dev.off()

# Get descriptives for Table 1
sum(dat$iv_changes_gov_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)
sum(dat$iv_changes_gains_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)
sum(dat$iv_changes_losses_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)

summary(abs(dat$iv_changes_gov)[which(dat$iv_changes_gov_sig==1)])
summary(abs(dat$iv_changes_gains)[which(dat$iv_changes_gains_sig==1)])
summary(abs(dat$iv_changes_losses)[which(dat$iv_changes_losses_sig==1)])

# Figure 3b
dat$Policy <- as.factor(dat$policy)
levels(dat$Policy) <- c("No","Yes")
dat$Policy <- relevel(dat$Policy,ref="Yes")
pdf("ivpublicsupport_imp.pdf")
ggplot(subset(dat,!is.na(dat$Policy)),aes(gw_fvr_dk,iv_changes_gov,colour=Policy)) + geom_point() + theme_bw() + geom_vline(xintercept=50,linetype="dashed") + geom_hline(yintercept=c(-1,0,1),linetype="dashed",colour=c("grey","black","grey")) + ylab("Total Issue Voting") + xlab("Public Support (%)") + scale_colour_grey(name="Policy Implemented") + theme(legend.position="bottom")
dev.off()

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==0]==1,na.rm=T)


# Make graph considering whether issue voting is significant or not (positive and negative). Figure 3a
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0)) + theme(legend.position="bottom")
dev.off()


# Table 2
mod1 <- lm(policy ~ iv_changes_gov_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*gw_fvr_dk  + totalchanges*gw_fvr_dk,data=dat)
library(apsrtable)
apsrtable(mod1,mod2,mod3,stars="default")

# Figure 4a
library(margins)
m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gov_sig"),]
pdf("absivme.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Figure 4b
m <- summary(margins(mod2,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("ogivme.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Gov IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()


# Figure 4c
m <- summary(margins(mod3,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_losses_sig"),]
pdf("goivme.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Gov-to-Opp IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()


# Try looking at interaction the other way around
m <- summary(margins(mod1,at=list(iv_changes_gov_sig=c(0,1))))
m <- m[which(m$factor=="gw_fvr_dk"),]
m$iv_changes_gov_sig <- as.factor(m$iv_changes_gov_sig)
levels(m$iv_changes_gov_sig) <- c("No","Yes")

# Figure A6a
pdf("absivme_rev.pdf")
ggplot(m,aes(iv_changes_gov_sig,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Issue Voting") + ylab("AME Public Support")
dev.off()

# Figure A6b
m <- summary(margins(mod2,at=list(iv_changes_gains_sig=c(0,1))))
m <- m[which(m$factor=="gw_fvr_dk"),]
m$iv_changes_gains_sig <- as.factor(m$iv_changes_gains_sig)
levels(m$iv_changes_gains_sig) <- c("No","Yes")
pdf("ogivme_rev.pdf")
ggplot(m,aes(iv_changes_gains_sig,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Issue Voting") + ylab("AME Public Support")
dev.off()

# Figure A6c
m <- summary(margins(mod3,at=list(iv_changes_losses_sig=c(0,1))))
m <- m[which(m$factor=="gw_fvr_dk"),]
m$iv_changes_losses_sig <- as.factor(m$iv_changes_losses_sig)
levels(m$iv_changes_losses_sig) <- c("No","Yes")
pdf("goivme_rev.pdf")
ggplot(m,aes(iv_changes_losses_sig,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Issue Voting") + ylab("AME Public Support")
dev.off()

# Now do analyses with high-income preferences 

# Individual data points. Figure A14
pdf("ivpublicsupport_highincome.pdf")
ggplot(dat,aes(pi90_fvr_dk,iv_changes_gov)) + geom_point() + theme_bw() + geom_hline(yintercept=0,linetype="dashed") + ylab("Total Issue Voting") + xlab("Public Support (%)")
dev.off()

# Figure A16b
dat$Policy <- as.factor(dat$policy)
levels(dat$Policy) <- c("No","Yes")
dat$Policy <- relevel(dat$Policy,ref="Yes")
pdf("ivpublicsupport_imp_highincome.pdf")
ggplot(subset(dat,!is.na(dat$Policy)),aes(gw_fvr_dk,iv_changes_gov,colour=Policy)) + geom_point() + theme_bw() + geom_vline(xintercept=50,linetype="dashed") + geom_hline(yintercept=c(-1,0,1),linetype="dashed",colour=c("grey","black","grey")) + ylab("Total Issue Voting") + xlab("Public Support (%)") + scale_colour_grey(name="Policy Implemented")
dev.off()

mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_gov_sig==1]==1,na.rm=T)
mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_gov_sig==0]==1,na.rm=T)

mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_gains_sig==1]==1,na.rm=T)
mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_gains_sig==0]==1,na.rm=T)

mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_losses_sig==1]==1,na.rm=T)
mean(dat$policy[dat$pi90_fvr_dk>50 & dat$iv_changes_losses_sig==0]==1,na.rm=T)

# 45 % opinion threshold
mean(dat$policy[dat$pi90_fvr_dk>45 & dat$iv_changes_gains_sig==1]==1,na.rm=T)
mean(dat$policy[dat$pi90_fvr_dk>45 & dat$iv_changes_gains_sig==0]==1,na.rm=T)


# Make graph considering whether issue voting is significant or not (positive and negative). Figure A16a
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$policy[which(dat$pi90_fvr_dk>i-20 & dat$pi90_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$policy[which(dat$pi90_fvr_dk>i-20 & dat$pi90_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$policy[which(dat$pi90_fvr_dk>i-20 & dat$pi90_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(pi90_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$pi90_fvr_dk_fvr_dk <- df$pi90_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_highincome.pdf")
ggplot(df,aes(pi90_fvr_dk_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Table A13
mod1 <- lm(policy ~ iv_changes_gov_sig*pi90_fvr_dk + totalchanges*pi90_fvr_dk,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*pi90_fvr_dk + totalchanges*pi90_fvr_dk,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*pi90_fvr_dk  + totalchanges*pi90_fvr_dk,data=dat)
library(apsrtable)
apsrtable(mod1,mod2,mod3,stars="default")

library(margins)
# Figure A17a
m <- summary(margins(mod1,at=list(pi90_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gov_sig"),]
pdf("absivme_highincome.pdf")
ggplot(m,aes(pi90_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Figure A17b
m <- summary(margins(mod2,at=list(pi90_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("ogivme_highincome.pdf")
ggplot(m,aes(pi90_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Figure A17c
m <- summary(margins(mod3,at=list(pi90_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_losses_sig"),]
pdf("goivme_highincome.pdf")
ggplot(m,aes(pi90_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()


# Control for Intensity
load("data/intensity.Rda")
dat <- merge(dat,dat2,by="org_var")

# Rescale intensity from 0 to 1
dat$intensity <- (dat$intensity-min(dat$intensity,na.rm=T))/(max(dat$intensity,na.rm=T)-min(dat$intensity,na.rm=T))

mod1 <- lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk + intensity*gw_fvr_dk,data=dat)


# Marginal efects of opp-to-govt IV controlling for intensity
m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("gainsimp_intensity.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Marginal effects of intensity
m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="intensity"),]
pdf("intensity.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Intensity") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()



# Run models with continuing IV variables
dat$iv_changes_gains[dat$iv_changes_gains<0] <- NA
dat$iv_changes_losses[dat$iv_changes_losses>0] <- NA

# Create dataset with gains sorted from largest to smallest
df <- data.frame(iv_changes_gains=dat$iv_changes_gains,issue=dat$issueEnglish,year=dat$year)
# Remove missing
df <- df[which(!is.na(dat$iv_changes_gains)),]
df <- df[order(df$iv_changes_gains,decreasing=T),]

# Table A6
mod1 <- lm(policy ~ iv_changes_gov*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod2 <- lm(policy ~ iv_changes_gains*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod3 <- lm(policy ~ iv_changes_losses*gw_fvr_dk  + totalchanges*gw_fvr_dk,data=dat)
apsrtable(mod1,mod2,mod3,stars="default")

m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gov"),]
pdf("absivme_cont.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

m <- summary(margins(mod2,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains"),]
pdf("ogivme_cont.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

m <- summary(margins(mod3,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_losses"),]
pdf("goivme_cont.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Govt-to-Opp IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Now run models with support among government party supporters. Table A8
mod1 <- lm(policy ~ iv_changes_gov_sig*govsupport + totalchanges*govsupport,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*govsupport + totalchanges*govsupport,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*govsupport + totalchanges*govsupport,data=dat)
apsrtable(mod1,mod2,mod3,stars="default")


# Total IV. Figure 5
m <- summary(margins(mod1,at=list(govsupport=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_gov_sig"),]
pdf("absgovtparty.pdf")
ggplot(m,aes(govsupport,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()

# Opp-to-Govt. Figure 5b
m <- summary(margins(mod2,at=list(govsupport=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("gainsgovtparty.pdf")
ggplot(m,aes(govsupport,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()

# Govt-to-Opp. Figure 5c
m <- summary(margins(mod3,at=list(govsupport=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_losses_sig"),]
pdf("lossesgovtparty.pdf")
ggplot(m,aes(govsupport,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Govt-to-Opp IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()



# Distributions of support for adopted policies among govt party supporters and overall population. Figure A7
df <- data.frame(support=c(dat$gw_fvr_dk,dat$govsupport),policy=c(dat$policy,dat$policy),group=c(rep("All Respondents",length(dat$gw_fvr_dk)),rep("Govt Party Supporters",length(dat$gw_fvr_dk))))
pdf("supportdistribs.pdf")
ggplot(df[df$policy==1,],aes(support,fill=group)) + geom_density(alpha=0.3) + theme_bw() + scale_fill_grey() + ylab("Density") + xlab("Support (%)") + theme(legend.title=element_blank())
dev.off()



# Post-election government party supporters
mean(dat$policy[which(dat$iv_changes_gov_sig==1 & dat$govsupport_post>50)],na.rm=T)
mean(dat$policy[which(dat$iv_changes_gov_sig==0 & dat$govsupport_post>50)],na.rm=T)

mean(dat$policy[which(dat$iv_changes_gains_sig==1 & dat$govsupport_post>50)],na.rm=T)
mean(dat$policy[which(dat$iv_changes_gains_sig==0 & dat$govsupport_post>50)],na.rm=T)

mean(dat$policy[which(dat$iv_changes_losses_sig==1 & dat$govsupport_post>50)],na.rm=T)
mean(dat$policy[which(dat$iv_changes_losses_sig==0 & dat$govsupport_post>50)],na.rm=T)


# Table A9
mod1 <- lm(policy ~ iv_changes_gov_sig*govsupport_post + totalchanges*govsupport_post,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*govsupport_post + totalchanges*govsupport_post,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*govsupport_post + totalchanges*govsupport_post,data=dat)
apsrtable(mod1,mod2,mod3,stars="default")

# Total IV. Figure A8a
m <- summary(margins(mod1,at=list(govsupport_post=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_gov_sig"),]
pdf("absgovtparty_post.pdf")
ggplot(m,aes(govsupport_post,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()

# Opp-to-govt. Figure A8b
m <- summary(margins(mod2,at=list(govsupport_post=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("gainsgovtparty_post.pdf")
ggplot(m,aes(govsupport_post,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Gov IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()

# Govt-to-opp. Figure A8c
m <- summary(margins(mod3,at=list(govsupport_post=c(10,20,30,40,50,60,70,80,90))))
m <- m[which(m$factor=="iv_changes_losses_sig"),]
pdf("lossesgovtparty_post.pdf")
ggplot(m,aes(govsupport_post,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Gov-to-Opp IV") + geom_rug(data=dat,mapping=aes(govsupport_post,0))
dev.off()

# Do make these separately for cases of new governments and incumbents
dat$newgov <- ifelse(dat$year %in% c(1976,1982,1991,1994,2006),1,0)   
dat$newgov <- as.factor(dat$newgov)
levels(dat$newgov) <- c("Incumbent","Non-incumbent")

# Create variable for significant gains/losses/not-significant
dat$sigchanges <- 0
dat$sigchanges[which(dat$iv_changes_gains_sig==1)] <- 1
dat$sigchanges[which(dat$iv_changes_losses_sig==1)] <- -1
dat$Significant <- as.factor(dat$sigchanges)
levels(dat$Significant) <- c("Govt-to-Opp","No Significant Change","Opp-to-Govt")
prop.test(table(dat$sigchanges[dat$newgov=="Non-incumbent" & dat$gw_fvr_dk>50],dat$policy[dat$newgov=="Non-incumbent" & dat$gw_fvr_dk>50])[2:3,2:1],alternative="less")
# Figure A9
pdf("ivpublicsupport_incumbent.pdf")
ggplot(dat,aes(gw_fvr_dk,iv_changes_gov,colour=Significant)) + geom_point() + theme_bw() + facet_wrap(~newgov) + ylab("Issue Voting") + xlab("Public Support (%)") + scale_colour_grey(start=0.8,end=0.2)
dev.off()



# Run jacknife
# drop missing
dat <- dat[!is.na(dat$iv_changes_gov),]
pvalues <- vector()
for(i in 1:nrow(dat)){
  mod1 <-lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat[-i,])
  pvalues <- c(pvalues,summary(mod1)$coefficients[5,4])
}

library(ggplot2)
df <- data.frame(number=1:202,pvalues=pvalues)
#Figure A11
pdf("pvalues_jackknife.pdf")
ggplot(df,aes(number,pvalues)) + geom_point() + coord_flip() + ylab("P-value") + xlab("") + theme_bw() + theme(axis.title.y=element_blank(),                                                                                                               axis.text.y=element_blank(),                                                                                                            axis.ticks.y=element_blank())
dev.off()


# Get issue numbers with p-values greater than 0.05
higher <- which(pvalues>=0.05)
plots <- list()
library(margins)
#start counter variable
j <- 1
for(i in higher){
  mod1 <-lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat[-i,])
  m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
  m <- m[which(m$factor=="iv_changes_gains_sig"),]
  p <- ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Gains") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
  plots[[j]] <- p
  j <- j+1 # increment counter variable
}


# Now get DFBETAs
mod1 <-lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
dfb <- as.data.frame(dfbeta(mod1))
dfbdf <- data.frame(year=dat$year,issueEnglish=dat$issueEnglish,gw_fvr_dk=dat$gw_fvr_dk,iv_changes_gains_sig=dat$iv_changes_gains_sig,policy=dat$policy)[which(rownames(dat) %in% rownames(dfb)),]
dfbdf$dfb <- dfb$"iv_changes_gains_sig:gw_fvr_dk"
dfbdf <- dfbdf[order(dfbdf$dfb,decreasing=T),]
# Are any above the cut-off?
sum(dfbdf$dfb>(2/sqrt(nrow(dat))))

df <- data.frame(number=1:151,dfbeta=dfb$"iv_changes_gains_sig:gw_fvr_dk")
#Figure A12
pdf("dfbetas.pdf")
ggplot(df,aes(number,dfbeta)) + geom_point() + ylab("DFBETA") + xlab("") + theme_bw() + theme(axis.title.x=element_blank(),
                                                                                              axis.text.x = element_blank(),
                                                                                              axis.ticks.x=element_blank()) + geom_hline(mapping=aes(yintercept=c(2/sqrt(151))),linetype="dashed") + geom_hline(mapping=aes(yintercept=c(-2/sqrt(151))),linetype="dashed")
dev.off()

# Identify points with high issue voting and public support and which were implemented
torm <- which(dat$gw_fvr_dk>60 & dat$iv_changes_gov>4)
# Now run the model without them. Table A10
mod1 <-lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat[-torm,])

library(apsrtable)
apsrtable(mod1,stars="default")

m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
# Table A13
pdf("ogivme_infl.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Gov IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()



# Social Democratic Governments
dat$socdem <- ifelse(dat$year %in% c(1956,1958,1960,1964,1968,1970,1973,1982,1985,1988,1994,1998,2002,2014),1,0)


mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==1 & dat$socdem==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==0 & dat$socdem==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==1 & dat$socdem==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==0 & dat$socdem==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==1 & dat$socdem==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==0 & dat$socdem==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==1 & dat$socdem==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==0 & dat$socdem==1]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==1 & dat$socdem==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==0 & dat$socdem==1]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==1 & dat$socdem==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==0 & dat$socdem==1]==1,na.rm=T)

# Minority governments
dat$notminority <- ifelse(dat$year %in% c(1956,1968,1976,1979,2006),1,0)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==1 & dat$notminority==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==0 & dat$notminority==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==1 & dat$notminority==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gains_sig==0 & dat$notminority==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==1 & dat$notminority==0]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==0 & dat$notminority==0]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==1 & dat$notminority==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_gov_sig==0 & dat$notminority==1]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>40 & dat$iv_changes_gains_sig==1 & dat$notminority==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>40 & dat$iv_changes_gains_sig==0 & dat$notminority==1]==1,na.rm=T)

mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==1 & dat$notminority==1]==1,na.rm=T)
mean(dat$policy[dat$gw_fvr_dk>50 & dat$iv_changes_losses_sig==0 & dat$notminority==1]==1,na.rm=T)

# Table A11
mod1 <- lm(policy ~ iv_changes_gov_sig*gw_fvr_dk + socdem*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + socdem*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*gw_fvr_dk + socdem*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
library(apsrtable)
apsrtable(mod1,mod2,mod3,stars="default")

# Table A12
mod1 <- lm(policy ~ iv_changes_gov_sig*gw_fvr_dk + notminority*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod2 <- lm(policy ~ iv_changes_gains_sig*gw_fvr_dk + notminority*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
mod3 <- lm(policy ~ iv_changes_losses_sig*gw_fvr_dk + notminority*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat)
library(apsrtable)
apsrtable(mod1,mod2,mod3,stars="default")



# Graphs with implementation by year

# Figure A10
# Make graph considering whether issue voting is significant or not (positive and negative)
# Year 0 
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$year0policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$year0policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$year0policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_year0.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Year 1 
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$year1policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$year1policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$year1policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_year1.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Year 2
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$year2policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$year2policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$year2policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_year2.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Year 3
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_year3.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Year 4
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$year3policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_year4.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())  + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

